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(57) Abstract: Systems and methods for solving Finite element models, wherein the matrix that governs the solution is modified by 
£J adjusting the weighting coefficients of the matrix so that the elements which lie on the diagonal of the matrix are non-negative and 

the elements which are off the diagonal are non-positive. In one embodiment, a system is discretized on a finite element mesh with 
Q the contribution of each node to the discretization being weighted based upon the direction of fluid flow across each element The 
^ nodes which are upstream from the other nodes of the respective elements are weighted more heavily to cause the resulting matrix 
^ to be substantially diagonal. This matrix is solved using traditional techniques. 
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METHOD AND SYSTEM FOR SOLVING FINITE ELEMENT MODELS 
USING MULTI - PHASE PHYSICS 

Related Applications 

[000 11 This application claims priority under 35 U.S.C. 
§ 119(e) to United States Patent Application 
No. 60/215,697 entitled "Method and System for Oil 
Reservoir Simulation and Modeling" by Stephen R. Kennon, 
Kok Thye Lim, Scott A. Canaan, Steven B. Ward, Stuart W. 
Pond, Jr. and Edward J. Barragy, filed June 29, 2000, 
which is incorporated by reference as if set forth in 
its entirety herein. 

Field of the Invention. 

[0002] The invention relates generally to methods for 

modeling physical systems using finite element analysis 
and, more specifically, to methods for solving finite 
element models wherein the matrix that controls the 
solution for the model is configured to maintain the 
properties of monotonicity and linearity preservation. 

Background of the Invention. 

[0003] Physical systems can be modeled mathematically to 

simulate their behavior under certain conditions. There 
are a wide variety of means to model these systems, 
ranging from the very simplistic to the extremely 
complicated. One of the more sophisticated means to 
model physical systems is through the use of finite 
element analysis. As the name implies, finite element 
analysis involves the representation of individual, 
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finite elements of the physical system in a mathematical 
model and the solution of this model in the presence of 
a predetermined set of boundary conditions. Another 
comparable means to model physical systems is through 
the use of finite difference analysis. Finite 
difference analysis involves the modeling of individual 
points within a modeled space and computing the 
differences between these points. Finite difference 
analysis is often used for simulating the dynamic 
behavior of fluids. 

[0004] Traditional finite difference techniques and finite 
element techniques using streamlined, upwinding methods 
are typically used to simulate the production of oil in 
a reservoir. While each of these techniques has its own 
advantages, it also has its own disadvantages. 
Generally speaking, the finite difference techniques 
produce physically realistic results, but they are not 
very accurate. The finite element techniques, on the 
other hand, are more accurate, but they produce results . 
which are not physically realistic. 

[0005] As a result of their respective disadvantages, both 
the finite difference techniques and finite element 
techniques which are conventionally used normally 
require a great deal of computer resources. In the case 
of finite difference techniques, reasonable accuracy can 
be achieved, but this requires many more nodes then 
would be necessary in a finite element model. This 
increases the amount of memory and CPU time which are 
needed to compute an accurate solution. In the case of 
finite element techniques, comparable accuracy can be 
achieved with fewer nodes, but the solutions may not be 
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realistic. For example, concentrations may be more than 
100 percent, or permeabilities may be negative. 
Consequently, there may be problems for which the 
solution does not converge, or for which the solution 
may converge very slowly. These techniques are 
therefore less reliable and may require a large number 
of iterations before acceptable accuracy is achieved. 

[0006] Because each of these standard techniques has its 
own drawbacks, and because these drawbacks increase in 
the amount of computer resources which are necessary to 
generate acceptable solutions, it would be desirable to 
provide a method for modeling systems such as oil 
reservoirs which reliably produces accurate, realistic 
solutions for these systems. 

Summary of the Invention 

[0007] One or more of the problems outlined above may be 
solved by the various embodiments of the- invention. 
Broadly speaking, the invention comprises systems and 
methods for solving finite element models, wherein the 
matrix that governs the solution is modified so that it 
satisfies the properties of monotonicity and 
preservation of linearity. This is accomplished by 
adjusting the weighting coefficients of the matrix so 
that the elements which lie on the diagonal of the 
matrix are non-negative and the elements which are off 
the diagonal are non-positive. This causes the solution 
to be both more accurate than traditional finite 
difference techniques and more realistic than 
traditional finite element techniques. 
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[0008] The present techniques may provide various 

advantages over the prior art. For instance, because 
the present techniques are more accurate, the number of 
iterations which are required to produce an acceptable 
result may be reduced, thereby conserving computing 
resources. The present techniques may also cause 
solutions to more reliably converge, since they are 
constrained to the range of realistic results and do not 
oscillate in the same manner as the solutions of 
traditional finite element techniques. The present 
techniques are therefore generally faster, more 
accurate, less resource -intensive and more reliable than 
prior art techniques. 

[0009] One embodiment of the present invention comprises a 
method for generating a solution for a modeled system 
using finite element techniques, wherein the system is 
discretized on a finite element mesh and wherein the 
contribution of each node to the discretization is 
weighted based upon the direction of fluid flow across 
each element. The contributions are weighted to favor 
the nodes which are upstream from the other nodes of the 
respective elements. Weighting the contributions of the 
nodes in this manner causes the resulting matrix to be 
essentially diagonal. That is, the elements which lie 
on the diagonal of the matrix are emphasized, while the 
elements which are off the diagonal are de -emphasized. 
Once the matrix has been modified in this manner, it is 
solved using traditional techniques. 



[0010] An alternate embodiment of the present invention 

comprises a method for modeling the flow of multiphase 
fluids in an oil reservoir. In this embodiment, a four 
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dimensional finite element mesh corresponding to the oil 
reservoir is constructed. A solution is generated for 
the mesh by distributing the residuals of a finite 
element operator among the nodes of the mesh and then 
assembling and solving the system using traditional 
finite element techniques. In this embodiment, 
distributing the residual of the finite element operator 
comprises, for each element, computing the value of the 
finite element operator for each node, dividing these 
values into a set of positive values and a set of 
negative values, distributing the residual to the nodes 
to emphasize the diagonal elements and de-emphasize the 
off -diagonal elements and then computing a tangent 
operator for the matrix. 

[0011] Another alternate embodiment of the present 

invention comprises a software product which is embodied 
in a computer readable medium. The software product 
comprises a plurality of instructions which are 
configured to cause a computer or any other data 
processor to perform the steps of a method in accordance 
with the present disclosure. The computer readable 
medium may comprise any of a number of different media, 
including but not limited to CD ROMs, DVDs, floppy 
disks, computer memories, magnetic tapes, etc. Still 
another alternate embodiment of the present invention . 
comprises a computer system which is configured to 
perform the steps of a method in accordance with the 
present disclosure. The computer system may be 
programmed to perform such a method by instructions 
which are contained in hardware, firmware or software. 
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[0012] Numerous additional alternative embodiments are also 
possible. 

Brief Description of the Drawings 

[0013] Other objects and advantages of the invention may 
become apparent upon reading the following detailed 
description and upon reference to the accompanying 
drawings . 

[0014] FIGURE 1 is a diagram illustrating simplices in two 
and three dimensions. 

[0015] FIGURE 2 is a diagram illustrating the difference 

between the actual values in the physical system and the 
solutions generated using traditional finite element 
techniques and techniques according to the present 
invention. 

[0016] FIGURE 3 is a diagram illustrating the flow of fluid 
across a first two-dimensional simplex element. 

[0017] FIGURE 4 is a diagram illustrating the flow of fluid 
across the first two-dimensional simplex element in a 
direction opposite the flow illustrated in FIGURE 3. 

[0018] FIGURE 5 is a diagram illustrating the flow of fluid 
across a second two-dimensional simplex element . 

[0019] FIGURE 6 is a diagram illustrating the flow of fluid 
across the second two-dimensional simplex element in a 
direction opposite the flow illustrated in FIGURE 5. 

[0020] FIGURE 7 is a diagram illustrating a two-dimensional 
projection of a time-space domain. 
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[0021] FIGURE 8 is a flow diagram illustrating a method in 
accordance with one embodiment of the present invention. 

[0022] While the invention is subject to various 
modifications and alternative forms , specific 
embodiments thereof are shown by way of example in the 
drawings and the accompanying detailed description- It 
should be understood, however, that the drawings and 
detailed description are not intended to limit the 
invention to the particular embodiment which is 
described- This disclosure is instead intended to cover 
all modifications, equivalents and alternatives falling 
within the scope of the present invention as defined by 
the appended claims. 

Detailed Description of a Preferred Embodiment 

[0023] A preferred embodiment of the invention is described 
below. It should be noted that this and any other 
embodiments described below are exemplary and are 
intended to be illustrative of the invention rather than 
limiting. 

[0024] Broadly speaking, the invention comprises methods 

for solving finite element models, wherein the solution 
is guaranteed both to be monotonic and to satisfy 
linearity preservation. This ensures that the solution 
will be both accurate and realistic. The mono tonicity 
and the linearity of the solution are ensured by 
generating a matrix for the solution of the finite 
element model such that the elements of the. matrix which 
are on the diagonal are non-negative, and the elements 
which are off -diagonal are non-positive. This matrix 
can then be solved using conventional finite element 
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techniques to generate the desired accurate, realistic 
solution. 

[0025] The present methods provide a number of advantages 
over traditional finite element and finite difference 
techniques. As indicated above, the present methods 
provide the accuracy of traditional finite element 
techniques without the disadvantages of unrealistic, 
oscillating solutions. Similarly, the present methods 
provide the stability of traditional finite difference 
techniques while providing much greater accuracy. 

[0026] The accuracy and realism provided by the present 

methods lead to additional advantages in the utilization 
of computer resources. For example, the present methods 
may allow for faster generation of solutions for 
physical systems . This is a result of several things . 
First, because the solution is physically realistic, it 
does not oscillate, and consequently tends to converge 
more rapidly. Further, because generation of the 
solution using the present methods typically requires 
less iterations and is faster than traditional finite 
element techniques, less CPU time is required to 
generate the solution. The present methods may also use 
less memory than traditional techniques. Because the 
present methods are generally more accurate than 
traditional finite difference techniques, is not 
necessary to use as many nodes to model the physical 
system. Consequently, less memory is required to store 
the model (and less CPU time is required to compute a 
solution based on these nodes.) the present methods may 
also be more reliable than traditional techniques. The 
oscillation of solutions generated using traditional 
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finite element techniques may, in some instances, 
prevent a solution from converging at all. As noted 
above, the present methods provide solutions which are 
realistic (e.g., they do not oscillate between realistic 
solutions and solutions which are not possible,) so the 
solutions will reliably converge. 

[0027] In finite element modeling, the region that is to be 
analyzed is broken up into sub-regions called elements. 
This process of dividing the region into sub-regions may 
be referred to as discretisation or mesh generation. 
The region is represented by functions defined over each 
element . This generates a number of local functions 
that are much simpler than those which would be required 
to represent the entire region. The next step is to 
analyze the response for each element. This is 
accomplished by building a matrix that defines the 
properties of the various elements within the region and 
a vector that defines the residuals acting on each 
element in the domain. Once all the element matrices 
and vectors have been created, they are combined into a 
matrix equation. After applying boundary conditions, the 
matrix equation can be solved to obtain unknown nodal 
responses. Intra-element responses can be interpolated 
from nodal values using the functions which were defined 
over each element. 

[0028] The following description of one embodiment of the 
present method relates to the modeling of oil 
reservoirs. This embodiment is intended to be 
exemplary, and it should be noted that the present 
methods may be applied to a variety of physical systems, 
such as flow over aircraft, weather prediction, flow 
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through internal combustion engines, and other 
industrial fluid flows. This list of potential 
applications of the present methods is intended to be 
illustrative rather than limiting. 

[0029] The details of a preferred embodiment will be set 
forth below. It may be helpful, however, to first 
define a few terms. 

[0030] A node is a point in space. In finite element 

modeling, nodes form the vertices of the elements which 
are modeled. The nodes also form part of a mesh of 
nodes and edges which define the boundaries between 
elements in the modeled space. 

[0031] An edge is a line between two nodes which form 

vertices of an element. The edges form part of the mesh 
which defines the boundaries between elements in the 
modeled space. 

[0032] A simplex is a spatial configuration of n dimensions 
determined by n + 1 points in a space of dimension equal 
to or greater than n. In other words, a— simplex is a 
geometric spatial element having the minimum number of 
boundary points necessary to enclose a space in a given 
number of dimensions. For example, in two dimensions, 
a simplex comprises a triangle, together with the 
interior area bounded by the triangle (see FIGURE 1.) 
Two points are insufficient to form a simplex in two- 
dimensional space because no area is bounded by the 
points (and the lines which interconnect them.) While 
four points may be sufficient to bound a two-dimensional 
area, they do not comprise the minimum number of 
boundary points by which the two-dimensional area can be 
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bounded. In three dimensions, a simplex comprises a 
tetrahedron, which is bounded by four vertices (see 
FIGURE 1.) In four dimensions, a simplex comprises a 
hypertetrahedron (sometimes referred to as a hypertet) 
having five vertices. 

[0033] A mesh is a collection of elements that fill a 

space. These elements are representative of a system 
which resides in that space. Because each element can 
be defined by a plurality of nodes and/or the edges 
between those nodes, a mesh may alternatively be 
considered a collection of nodes and/or the edges 
between them. At various points in this disclosure, 
"mesh" will be used to alternately refer to collections 
of elements or nodes/edges, depending upon the context 
in which the term is used. The mesh may also be 
referred to herein as a finite element model or simply a 
model . 

[0034] There is a fundamental difference between modeling 

systems in which there is a single fluid phase and those 
in which there are multiple phases. When there are 
multiple fluid phases, there is a convection and 
potentially a mixing of the phases. (It is assumed that 
there is a pressure which drives the fluids and causes 
the front between them to be transported.) The finite 
element model must therefore include convective 
operators or components. These operators would not be 
necessary if only one fluid phase were modeled, or if 
the multiple fluid phases were completely mixed or 
emulsified. 
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[0035] In a system having multiple fluid phases, 

characteristics of the system may change very abruptly. 
This situation is illustrated in FIGURE 2. FIGURE 2 is 
a diagram illustrating the difference between the actual 
values in the physical system and the solutions 
generated using traditional finite element techniques 
and techniques according to the present invention. It 
can be seen that the actual ("real world") values 11 
form what is essentially a step function, changing from 
a first value vl to a second value vO immediately upon 
reaching a distance dl. In this example, the values vO 
and vl define the limits of the range of physically 
realizable values. 

[0036] When traditional finite element techniques are used 
to model such a system, the solution will typically tend 
to oscillate. It can be seen in FIGURE 2 that the 
solution 12 generated using these techniques oscillates 
above and below the actual physical values which it is 
attempting to model. Two of the most significant 
problems with this solution are the obvious inaccuracy 
(the solution does not closely track the actual values) 
and the fact that, at some points, the solution is 
greater than the maximum physically realizable value or 
less than the minimum physically realizable value. The 
oscillations which occur in solutions generated using 
traditional finite element techniques may therefore 
produce unrealistic results, or may even fail to 
converge at all . 

[0037] Methods in accordance with this disclosure, however, 
do not suffer from these disadvantages because the 
solutions are constrained to be realistic. This 
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constraint is achieved by ensuring that the matrix which 
governs the solution vector (i.e., the Jacobian matrix) 
satisfies the properties of monotonicity and 
preservation of linearity. In more practical terms, 
this consists of ensuring that all of the coefficients 
on the diagonal of the matrix are non-negative and that 
all of the coefficients off the diagonal of the matrix 
are non-positive. When the matrix has these 
characteristics, the solution remains realistic at all 
points, and more accurately tracks the behavior of the 
physical system, as indicated by curve 13 in FIGURE 2. 
The model is therefore more likely to produce accurate 
results and to reliably converge. 

[003 8] The Jacobian matrix is derived from the partial 
differential equations that describe multiphase flow 
through porous media based on Darcy's law, and 
representing conservation of mass and momentum in the 
system. These equations are discretized onto the finite 
element mesh, and the coefficients of the equations form 
a matrix. Each coefficient will be either positive or 
negative, depending upon whether the gradient of the 
equation is positive or negative (i.e., whether the 
fluid flow is positive or negative) . If the finite 
element technique is applied to the differential 
equations, a large set of algebraic equations is 
produced. That set of algebraic equations is nonlinear. 
The set of equations can be solved by using Newton's 
method (which involves taking the first derivative of 
the equations) . The first derivative is the Jacobian 
matrix. The coefficients of this matrix are weighted to 
obtain a matrix which has non-negative terms on its 
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diagonal and non-positive terms off the diagonal. The 
appropriate weighting of the terms ensures the monotonic 
and linear properties of the matrix. 

[0039] The coefficients of the Jacobian matrix are 

essentially weighting factors for each of the nodes of 
the elements in the finite element mesh. In traditional 
finite element techniques, the nodes of an element are 
evenly weighted. For example, in a two-dimensional 
element, each of the three nodes would be weighted by 
1/3. If a three-dimensional element were used, each of 
the four nodes of the element would be weighted by 1/4. 
In the present methods, however, the nodes are weighted 
to reflect the flow of the fluids in the modeled system. 
This is illustrated with reference to FIGURES 3-6. 

[0040] FIGURE 3 shows a single, two-dimensional element (a 
triangle) . The flow of fluid with respect to the 
element is indicated by the arrows, which show that the 
fluid is flowing from the upper, right-hand portion of 
the figure to the lower, left-hand portion of the 
figure. Using traditional finite element techniques, 
each of nodes 21-23 would be evenly weighted (i.e., 
1/3, 1/3, 1/3). Using the present techniques, the nodes 
are not evenly weighted- Because the fluid is flowing 
generally from nodes 22 and 23 to node 21, each of nodes 
• 22 and 23 is weighted by 1/2, while node 21 is weighted 
by 0. Referring to FIGURE 4, it can be seen that the 
weighting of the nodes is different when the direction 
of the fluid flow changes. In FIGURE 4, the fluid flow 
is generally from the lower, left-hand portion of figure 
to the upper, right-hand portion of the figure. 
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Consequently, node 21 is weighted by 1, while nodes 22 
and 23 are weighted by 0. 

[0041] Referring to FIGURES 5 and 6, it can be seen that 

the weighting of the nodes is also affected by the shape 
of the corresponding element. In FIGURE 5, the flow of 
fluid is approximately perpendicular to an edge of the 
element. This situation is similar to that of FIGURE 3, 
in which each of the nodes which terminates the forward 
edge of the element was weighted by 1/2. In that case, 
the remaining node of the element was downstream from a 
midway point on the forward edge of the element. In 
FIGURE 5, the remaining node is instead directly 
downstream from one of the first to nodes . In this 
situation, the node 32 which is upstream from the last 
node 33 is weighted by 2/3, while the other upstream 
node 31is weighted by 1/3. In FIGURE 6, there is a 
single upstream node 33. Nodes 31 and 32 are both on 
the trailing edge of the element, so they are both 
weighted by 0. The upstream node 33, the other hand, is 
weighted by 1 . 

[0042] In each of these cases, the weighting of the nodes 
with respect to the direction of fluid flow effectively 
emphasizes the diagonal elements of the matrix, while 
reducing or eliminating the effect of the off-diagonal 
elements. This, in turn, ensures that the matrix has 
the properties of monotonicity and linearity. 

[0043] Although the foregoing is sufficient to enable a 

person of ordinary skill in the field of finite element 
analysis to implement the present techniques, a brief 
technical discussion of the techniques follows. 
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[0044] Considering a physical system which is to be modeled 
over time, a finite element mesh is first constructed 
(i.e., the system is discretized in the space-time 
domain) . Referring to FIGURE 7, a two-dimensional 
projection of the time-space domain is illustrated. 

[0045] ^ = ^space x fO» T ] i s the space-time domain to be 

modeled, where 

[0046] & is broken down (segmented) into a plurality of 
non-overlapping sub-domains . Each of these sub- 



u Q e = Q ee {1,...,N B } 

domains and « comprises a four- 

dimensional simplex (a hypertetrahedron) . A linear 
finite element basis is formed on each element (^e) , and 
trial functions and test functions are defined on each 
element . 



[0047] The flow of multi-phase fluids in porous media is 
governed by the equation: 
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+ V«(V i X i ) = 0 



Where: 



i = 1,2, ...,nphases 
V i --K(VP i + p i gVZ) 
P, -P 

Pi = P-P c ,i=2, ...,nphases 

Pi = Viscosity of i'th phase 

Bj = Formation Volume Factor of i'th phase 
Pi = Density of i'th phase 

<|> = Rock Porosity 
Z = Negative of Depth 
= Saturation of i'th phase 



Z Si"* 1 



[0048] 



[0049] 



A standard Galerkin Finite Element discretization of 




this multi-phase flow equation gives, for each element: 



[0050] 



[0051] 



The basic method of employing the present techniques 
is defined by the way that h is defined. The present 
method is therefore defined as follows: 



WO 02/03103 



18 



PCT/US01/20742 



f / - (j^ + i o /«i&«)(^) k *<«j'w. 

[0052] 

[0053] Referring to FIGURE 8 # a flow diagram illustrating a 
method in accordance with one embodiment of the present 
invention is shown. As indicated above, when the 
differential equations are discretized on the finite 
elements , they are basically converted from differential 
equations to algebraic relationships between the values 
of the nodes. This is accomplished by assuming that 
there is a linear variation of the corresponding 
property (e.g. saturation, pressure, permeability, 
prostate, etc.) between the nodes. This restricts the 
space of all possible solutions to only those solutions 
which are linear. When this assumption is used, a set 
of relations are derived from the differential equations 
for values of the coefficients which are linear over 
each element. When this is done, each element 
contributes a certain amount of the discretization to 
each one of its nodes. The contributions which are 
summed at each node are the residuals. It is desired to 
make the residuals equal 0 in order to solve the 
nonlinear algebraic system. 

[0054] To distribute the residuals of the finite element 

operator ( J } , Kj is first computed for each node of the 
4-tet (hypertetrahedron) . The resulting Kj are divided 
into positive (Kj + ) and negative (Kj") values. The 
T- 

residual J is then distributed to the nodes, and the 
tangent operator (Tjk) is computed. The system is then 
assembled and solved in the same manner as if an 
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ordinary finite element technique had been used. (Since 
the process by which finite element systems are 
assembled and solved is well known, it will not be 
described in further detail here) . 

[0055] The benefits and advantages which may be provided by 
the present invention have been described above with 
regard to specific embodiments. These benefits and 
advantages, and any elements or limitations that may 
cause them- to occur or to become more pronounced are not 
to be construed as a critical, required, or essential 
features of any or all of the claims. As used herein, 
the terms "comprises," "comprising," or any other 
variations thereof, are intended to be interpreted as 
non- exclusively including the elements or limitations 
which follow those terms. Accordingly, a process, 
method, article, or apparatus that comprises a list of 
elements does not include only those elements but may 
include other elements not expressly listed or , inherent 
to the claimed process, method, article, or apparatus. 

[0056] While the present invention has been described with 
reference to particular embodiments, it should be 
understood that the embodiments are illustrative and 
that the scope of the invention is not limited to these 
embodiments. Many variations, modifications, additions 
and improvements to the embodiments described above are 
possible. Particularly, these variations may include 
computers or other data processing devices, computer 
readable media (such as floppy disks, CD-ROMs, DVD-ROMs, 
etc.,) storage devices, computer memories and the like 
which contain software, firmware or other programming 
embodying the foregoing methods. It is contemplated 
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that these variations, modifications, additions and 
improvements fall within the scope of the invention as 
detailed within the following claims. 
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What is claimed is : 

1 . A method for solving a finite element model corresponding 
to a system in which there is a mult i -phase fluid flow 
comprising: 

generating a finite element matrix corresponding to the 
model, wherein the matrix contains a plurality of 
coefficients ; 

adjusting the coefficients to obtain a matrix in which 
on-diagonal elements are non-negative and off- 
diagonal elements are non-positive; and 

generating a solution for the model using the matrix 
using finite element techniques. 

2. The method of claim 1 wherein adjusting the coefficients 
comprises weighting the nodes of each element according to a 
direction of fluid flow across the element . 

3 . The method of claim 2 wherein weighting the nodes of each 
element according to a direction of fluid flow across the 
element comprises determining the direction of fluid flow 
across the element and weighting each node more heavily if the 
node is upstream from the other nodes of the element and less 
heavily if the node is downstream from the other nodes of the 
element . 

4. The method of claim 3 wherein each node is weighted 
more heavily if a greater portion of the element is downstream 
from the node than from other nodes of the element and less 
heavily if a smaller portion of the element is downstream from 
the node than from other nodes of the element . 
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5 . The method of claim 1 wherein the finite element matrix 

corresponds to a system in which there are at least two fluid 
phases. 

6. The method of claim 1 wherein the finite element matrix 
corresponds to a system in which there are three or more fluid 
phases . 

7 . The method of claim 1 wherein the finite element matrix 
corresponds to a four -dimensional finite element model. 

8 . The method of claim 1 wherein the system corresponds to 
an oil reservoir. 

9 . The method of claim 1 wherein the matrix is configured to 
produce a solution which is not physically unrealistic at any 
time. 

10. The method of claim 1 wherein the matrix is configured to 
produce a solution which is non-oscillating. 

11. The method of claim 1 further comprising discretizing a 
model of the system to produce a finite element mesh and 
generating the matrix based on the finite element mesh. 

12. A method for obtaining improved accuracy in solving 
finite element models comprising: 

discretizing a model of a system in which there is a 
multi-phase fluid flow; 
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generating a finite element matrix corresponding to the 

model, wherein the matrix is configured to produce a 
solution which is monotonic and preserves linearity; 
and 

generating a solution for the model using the matrix. 

13. A computer-readable medium which contains a plurality of 
instructions, wherein the instructons are configured to cause 
a computer to perform the method for solving a finite element 
model corresponding to a system in which there is a multi- 
phase fluid flow comprising: 

generating a finite element matrix corresponding to the 
model, wherein the matrix contains a plurality of 
coefficients ; 

adjusting the coefficients to obtain a matrix in which 
on-diagonal elements are non-negative and off- 
diagonal elements are non-positive; and 

generating a solution for the model using the matrix 
using finite element techniques. 

14. The computer-readable medium of claim 13 wherein 
adjusting the coefficients comprises weighting the nodes of 
each element according to a direction of fluid flow across the 
element . 

15. The computer- readable medium of claim 14 wherein 
weighting the nodes of each element according to a direction 
of fluid flow across the element comprises determining the 
direction of fluid flow across the element and weighting each 
node more heavily if the node is upstream from the other nodes 
of the element and less heavily if the node is downstream from 
the other nodes of the element. 
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16. The computer- readable medium of claim 15 wherein each 
node is weighted more heavily if a greater portion of the 
element is downstream from the node than from other nodes of 
the element and less heavily if a smaller portion of the 
element is downstream from the node than from other nodes of 
the element . 

17. The computer- readable medium of claim 13 wherein the 
finite element matrix corresponds to a system in which there 
are at least two fluid phases . 

18. The computer-readable medium of claim 13 wherein the 
finite element matrix corresponds to a system in which there 
are three or more fluid phases. 

19. The computer-readable medium of claim 13 wherein the 
finite element matrix corresponds to a four -dimensional finite 
element model . 

20. The computer- readable medium of claim 13 wherein the 
system corresponds to an oil reservoir. 

21. The computer-readable medium of claim 13 wherein the 
matrix is configured to produce a solution which is not 
physically unrealistic at any time. 



22 . The computer- readable medium of claim 13 wherein the 
matrix is configured to produce a solution which is non- 
oscillating. 
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23. The computer- readable medium of claim 13 wherein the 
method further comprises discretizing a model of the system to 
produce a finite element mesh and generating the matrix based 
on the finite element mesh. 
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Multi-Phase Flow 



Compute Kj for each 

node of the 4-Tet, 

j = 1,2, 3, 4, 5 Where: 

^--JVNj.K^Pt + PigVZ))^ 
♦ j ~ 1> • 5 



T 



Divide Kj into positive and nega- 
tive values 

K.J - KjifKjSO.elseO 
KJ = KjifKj^O.elseO 



T 



Distribute residual to the nodes 
according to 



E 1 ^ 



J 



Algorithm 1. To distribute residuals of 
the finite element operator C^j^spaiiai 



Compute tangent operator ■ 



Assemble system and solve. 
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